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REGENERATIVE  SIMULATION  WITH  INTERNAL  CONTROLS 

by 
Donald  L.  Iglehart  and  Peter  A.  W.  Lewis 

1 .   Introduction  and  Summary 

Simulators  are  frequently  faced  with  the  task  of  esti- 
mating a  parameter  associated  with  the  limiting  distribution  of 
a  stable  stochastic  process  which  is  being  simulated.   A 

methodology  based  on  regenerative  processes  for  obtaining  point 
estimates  and  confidence  intervals  for  such  parameters  from  a 
single  realization  of  the  process  has  recently  been  developed 
in  Crane  and  Iglehart  (1974a, b),  (1975a, b)  and  Iglehart  (1975), 
(1976a, b) ,  and  (1977) .   In  this  paper  we  shall  introduce  a  new 
technique,  internal  control  variables,  which  can  be  used  in  con- 
junction with  the  regenerative  method  for  obtaining  additional 
variance  reduction  for  the  estimates . 

Suppose  regenerative  process   {X  :t  >  0},  which  is  stable 
in  the  sense  that   X.  =»  X   as   t  ■>  °°   (here   =>   denotes  weak 
convergence) ,  is  being  simulated  by  generating  a  single  realization 
of  the  process.   For  convenience  think  of  this  process  as  being 
either  a  positive  recurrent  Markov  chain  or  the  waiting  time 
process  for  a  single  server  queue  with  traffic  intensity  less 
than  one.   Then  simulators  are  frequently  interested  in  estimating 
r  =  E{f(X)},  for  a  given  function   f.   The  principal  goal  of  the 
regenerative  method  is  to  produce  a  confidence  interval  estimate 
of   r.   The  regenerative  method  begins  by  observing  the  pairs 

1 


of  random  variable?   {  (Yi,/Tk)  :1  <_  k  <_  n}  ,  where   x,   is  the 
length  of  the  kth  regenerative  cycle  and   Y,   is  the  area  under 
the  function   f (X  )   in  the  kth  cycle.   Two  basic  facts  are  crucial 
for  the  regenerative  method  of  simulating  stable  processes.   First, 
the  pairs   {(Y,  ,x,)  :1  <_   k  <_  n}   are  independent  and  identically 
distributed  (i.i.d.),  and  second   r  =  E{f(X)}  =  EfY^/E^}. 
This  last  fact  suggests  that  a  natural,  strongly  consistent  point 

—  —  —In 

estimate  for   r   is   f(n)  =  Y(n)/x(n),  where   Y(n)  =  n    Ev-i  Yir 
and   x  (n)  =  n    lC--\  Tv  •   To  form  a  confidence  interval  for   r 
we  can  use  the  central  limit  theorem  (c.l.t.) 


(1.1)  /n  (r(n)  -  r)/(a/E{x1})  -  N(0,1)  , 


where   N(0,1)  denotes  a  mean  zero,  variance  one  normal  random  variable 

2  2 

and   a   =  E{ [Y,  -  rx,  ]  }.   Of  course,  the  constant   a/E{x,}   will 

have  to  be  estimated  in  most  simulations. 

The  idea  behind  internal  control  variables  is  to  introduce 
a  third  sequence  of  i.i.d.  random  variables   {C,:l  <_  k  <_   n}   which 
has  the  property  that   C,   is  defined  in  terms  of  the  kth  cycle 
and   E^cv^   is  known  (or  can  be  calculated  analytically) .   Then 
another  strongly  consistent  estimate  for   r   is 

n_1  S-i  [Yk +  3(ck  -  E<ck})] 

(1.2)  rCT(n)  E  ^ * 

x  (n) 

The  subscript  CT  is  meant  to  denote  "controlling  the  top"  in  the 

ratio  estimator;  similar  estimates  will  be  defined  for  "bottom  control." 


A  c.l.t.  analogous  to  (1.1)  also  holds  for   r     with   a   replaced 
by   a',  say.   Since  we  are  still  free  to  select   8,  we  choose  to 
do  so  in  such  a  way  as  to  minimize   o' .   Having  done  that  we  find 
that 

(a')2  =  a2[l  -  p2(C1/Y1  -  rTx)]  . 

To  obtain  significant  variance  reduction  a  control,  C,  ,  must  be 
found  which  is  highly  correlated  with  Y,  -  rx, ,  a  task  that  is 
not  always  easy. 

Section  2  of  the  paper  develops  the  method  of  internal  control 
variables  in  detail  and  contains  specific  examples  associated  with 
the  single  server  queue  and  Markov  chain  models.   These  ideas  are 
then  illustrated  with  simulation  results  and  numerical  calcula- 
tions in  Section  3.   Concluding  remarks  are  made  in  Section  4. 
In  particular  we  discuss  the  possibility  of  combining  internal 
controls  with  external  controls  to  obtain  further  variance  reduction. 


2 .   Internal  Controls 

2.1.   General  Ideas 

Let   X  =  {X  :t  ^  0 }  be  the  regenerative  process  being 

simulated.   Recall  that  a  process  is  regenerative  if  a  renewal 

j 
process   T  =  {T  :n  >_  0}   is  defined  on  the  same  probability  space 

as   X   and  if  the  portions  of  the   X   process  between  consecutive 

regeneration  points,  the   T  's,  are  i.i.d.   Typically  the   T  's 

n  n 

denote  the  times  the   X   process  enters  a  fixed  state  and  at 
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these  times  the  process  starts  over  from  scratch  independent  of 
its  past  history.   For  a  formal  definition  of  regenerative  process 
and  general  background  on  the  regenerative  method  of  simulation 
analysis  consult  Crane  and  Iglehart  (1975a) .   Let   TQ  =  0   for 
simplicity  and  define   tr  =  TR  -  T^,  k  >  1.   The  portion  of   X 
in  the  interval   [Tk_lfTk)   is  referred  to  as  the  kth  cycle  and 
is  of  length   x,.   Suppose  now  that   E{t1>  <  °°   and  that  the 
distribution  function  of   t1   is  aperiodic  (e.g.  its  support 
is  not  contained  in  a  set  of  the  form   {0 ,h, 2h, . . . } ,  where 
h  >  0)  .   Then  subject  to  mild  regularity  conditions   Xfc  =>  X 
as   t  /*  °°,  where   -  means   P{Xfc  <_   x>  ♦  P{X  <_  x}   for  all 
continuity  points  (x)  of  the  limit  distribution.   A  similar  result 
holds  when   t,   is  periodic;   see  Crane  and  Iglehart  (19  75a) .   The 
random  variable   X   is  frequently  thought  of  as  the  steady-state  con- 
figuration of  the  system  being  simulated.   Suppose   f   is  a 
measurable  function  from  the  state  space  of   X   to  the  real  line 
and  that  we  are  interested  in  estimating   r  =  E{f(X)}.   Define 
the  sequence  of  random  variables   (r.v.'s)   ^yt,:^  2l  ^   ky 

Tk 
(2.1)  Y,  =   /   f  (X  )ds  ,      k  >  1  . 

K       m  S  — 

ik-l 

If  the  time  parameter  of   X   is  discrete  (as  in  a  discrete  time 
Markov  chain),  the  integral  in  (2.1)  should  be  replaced  by  a  sum. 
The  regenerative  stucture  of   X,  the  process  being  simulated,  gives 
us  the  two  important  properties  stated  in  Section  1:   the  pairs 


{(Yk'Tk):1  -  k  -  n}   are  i-i-d-  and   r  =  E{Y1}/E{t1),  provided 
E  {  |  f  (X)  |  }  <  =o,  which  we  shall  always  assume.   The  c.l.t.  for  the 
ratio  estimator   r(n)  =  Y(n)/T(n)   indicated  in  (1.1)  is  proved 
in  Crane  and  Iglehart  (1975a)  and  follows  from  the  classical 


c.l.t.  for  the  i.i.d.  mean  zero,  finite  variance  r.v.'s 


2        ? 

Zk  =  Yk  -  rxk,  k  >_  1 .   We  always  assume   0  <  a   =  E{Z.~}  <  °° 


2  ^ 

The  variance  of  Z^,    o    ,  is  also  related  to  the  variance  of   r(n) 

through  the  asymptotic  relation  (as   n  ■+  °°) 


(2.2)  a2{r(n)}  -  n  *  (a2{  Z-^/E2^} )  . 


For  a  derivation  of  this  result  see  Cramer  (1964,   p.  354, 
eq.  (27.7.3)) .  A  number  of  point  estimates  and  confidence  intervals 
have  been  proposed  for  ratios  such  as   r   (see  Iglehart,  1975) ,  but 
the  simplest  conceptually  and  computationally  are  the  so-called 
classical  ones.   The  point  estimate  is   r(n)   and  the  confidence 
interval  is 


(2.3)  I(n)  =  [r(n)  -  z±_   /2s/~ '    r(n)  +  21-  /2s/t]     ' 

where   z,   ,^    =    <J>~  (1  -  y/2)  ,  the  inverse  of  the  standard  normal 
l-y/2 

distribution  function,  and   s   is  the  classical  point  estimate 

of   a   which  is  constructed  as  follows.   Let   s,,   [resp.   ^22] 

be  the  sample  variance  of  the   Y,  * s   [resp.  x,  ' s]   and   s12 

the  sample  covariance  of  the   Y,  ' s  and   Tk's-   Then   s   is  defined 

by 

S  =  [s11  -  2(Y/x)s12  +  (Y/T)2s22]1/^  . 
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Introduction  of  the  internal  control  variables  mentioned 

in  Section  1  is  done  with  the  hope  of  reducing  the  variance  of 

the  point  estimator  of   r.   This  in  turn  will  either  reduce  the 

number  of  cycles  that  need  to  be  simulated  for  fixed  precision 

or  reduce  the  length  of  the  confidence  interval, (2.3),  for  r. 

The  sequence  of  internal  control  variables   {C,  :k  >_  1}   is  defined 

on  the  same  probability  space  supporting  the  process   X   in  such 

a  way  that   C,   only  depends  on  X   (or  underlying  r.v.'s 

defining   X)  in  the  kth  cycle.   This  construction  of  the   C,  ' s 

will  insure  that  they  are  i.i.d.  because  of  the  basic  i.i.d. 

structure  of  the  cycles  of  a  regenerative  process.   The  control 

variables   C.   allowed,  however,  are  further  restricted  to  those 
k 

for  which  the  mean,  E{C,},  is  either  known  or  can  be  calculated 

analytically.   Thus  one  of  the  prices  we  must  be  prepared  to 

pay  to  obtain  variance  reduction  for  our  simulation  is  the 

analytical  work  of  computing   E{C, } .   Having  defined  the   C,  ' s 

we  proceed  to  form  new  ratio  estimators  for   r.   Starting  with 

the  ratio  estimator   r(n)  =  Y(n)/~(n)   we  can  either 

control  the   Y,  's,  the   x  's,  or  both  if  we  introduce  a  second 
K  k 

sequence  of  control  variables.   The  estimator   rp  (n) ,  proposed 
in  (1.2),  involves  controlling  the   Y,  ' s  only.   Here  the  sub- 
script CT   is  meant  to  suggest  "controlling  the  top."   For 
convenience  let  (for   k  >  1) 


and 


Yk  -  Yk  +  B(ck  -  Etck}) 


z  •  =  y;  -  n.  . 
k     k      k 


Observe  that  the   Y' ' s  are  i.i.d.  with   E{Yj }  =  E{Y  }  and  con- 
sequently the   ziVs  are  also  i.i.d.  with   E{Z'}  =  0.   Hence  the 
standard  c.l.t.  yields  (as   n  ■*■   °°) 

Y  z:/a{Z'}n1/2   =  n1/2[(Y/T)-r]/(a{Z1,}/:?)  -N(0,1) 
k=l   K     L  L 


which,  upon  replacing   t   by   E{t,}   in  the  denominator  (this 
can  be  justified  by  a  continuous  mapping  argument) ,  oecomes 

1/2  r*       i    \  i 

n         [rCT(n)    -   r] 

(c{Zj4/E{Tli)         "   N(0'1)' 

The  variance  of   Z ' ,  which  obviously  depends  on   3,  is  easily 
seen  to  be 


a2{Zp  =  a2{Z1}+26  covfZ-^C^}  +  ^a2^} 


2 
Now  select   3   so  as  to  minimize   a  {Z'K    This  yields 

cov{Z  ,C, } 
(2.4)  6*  =  - 


a2{C1) 


and 


(2.5)  a2{Z|}  =  a2{Z1)  [1  -  p2{Z1/C1}]  , 


where   p{Z,,C, }   is  the  correlation  coefficient  between   Z, 

and   C  .   If  we  use  the  c.l.t.  for  the  estimator   ^CT(n)  as  a  basis  for 


constructing  confidence  intervals  for   r,  then  we  want  a  control 

2 

variable   C,   which  will  minimize   a  {Z'}.   From  (2.5)  we  see 

2 

that  we  need  a  large  value  of   p  {Z,,C..}.   Since  the  length 

1  /2 
of  such  a  confidence  interval   is  proportional  to   o{Z'}/n  '     , 

to  reduce  the  number  of  cycles  required  by  a  factor  of  four  we 

1  /2 
need    |p{Z1/C1>|  >  (0.75)  7   =  0.866.   To  obtain  a  control,  C  , 

which  is  highly  correlated  with  either   Y,   or   t,   is  relatively 

easy  to  do.   However,  because   Y,   and   t,   are  themselves  highly 

correlated,  it  is  much  more  difficult  to  find  a  control   C. 

which  is  highly  correlated  with   Z, .   In  general,  we  shall  try 

to  find  a  control   C,   which  mimics   Z.   but  for  which  we  can 

still  calculate  its  mean,  e(c  ) .   As  usual  we  try  to  do  as  much 
analytically  as  possible,  before  having  to  resort  to  simulation. 
We  illustrate  some  possible  candidates  for  controls  in  the  context 
of  GI/G/1  queues  and  discrete  time  Markov  chains. 


2.2.   Internal  controls  for  GI/G/1  Queues 

In  the  GI/G/1  queue  we  assume  the  zeroth  customer  arrives 
at  time   t_  =  0,  finds  a  free  server,  and  experiences  a  service 

time   v_ .   The  nth  customer  arrives  at  time   t    and  experiences 

0  n        c 

a  service  time   v  .   Let  the  interarrival  times   t   -  t   .  =  u  , 

n  n    n-1     n 

n  >  1.   Assume  the  two  sequences   (v  :n  >  0}   and   {u  :n  >  1} 
—  n    —  n   — 

each  consists  of  independent,  identically  distributed  (i.i.d.) 
random  variables  (r.v.'s)  and  are  themselves  independent.   Let 

E{v  }  =  y  ,    E{u  }  =  A   ,  and  p  =  A/y,  where   0  <  A,y  <  °°.   Thus 


\i      has  the  interpretation  of  the  mean  service  rate  and   A   has 
the  interpretation  of  the  mean  interarrival  rate.   The  parameter 
p   is  called  the  traffic  intensity  and  is  the  natural  measure 
of  congestion  for  this  system.   We  shall  assume  that   p  <  1,  a 
necessary  and  sufficient  condition  for  the  system  to  be  stable. 
While  many  characteristics  of  interest  can  be  estimated 
using  the  regenerative  method,  we  shall  restrict  our  attention 

to  the  waiting  time  of  the  n —  customer,  W    (time  from  arrival 

3  n 

to  commencement  of  service) .   For  further  discussion  of  the 

simulation  of  the  GI/G/1  queue  see  Crane  and  Iglehart  (1974b) . 

To  obtain  a  representation  for  the  process   {W  :n  >^  0}   let 

X   =  v   ,  -  u   and  set   S~  =  0,  S   =X,  +»*»+X,n>l. 
n    n-1     n  0       n     1  n    — 

The  following  well-known  recursive  relationship  exists  for 

the   W  '  s: 

n 


Wn  =  0,    W  _,_,  =  [W   +  X  ,,]  +  ,    n  >  0 
0         n+1     n    n+1        — 


By  induction,  we  also  have 


W   =  max{S   -S.:k=0,  1,  ...,n},   n>0. 
n        n    k  — 


Using  the  strong  Markov  property  of  the  process   {S  :n  >_  0} 
it  can  be  shown  that  there  exists  a  sequence  of  integer-valued 

r.v.'s   {Tk:k  >  °)   sucn  tnat  T0  =  °'  Tk  <  Tk+1'  and  WT  =  ° 
with  probability  one.   In  other  words,  the  customers  numbered   T, 


are  those  lucky  fellows  who  arrive  to  find  a  free  server  and 
experience  no  waiting  in  the  queue.   The  fact  that  there  exists 
an  infinite  number  of  such  customers  in  the  GI/G/1  queue  is  a 
direct  consequence  of  the  assumption  that    p  <  1   and  the 
strong  law  of  large  numbers.   Thus   {W  :n  >^  0}   is  a  regenerative 
process.   If  we  let    t.  =  T   -  T,  _,  ,  k  >_  1 ,  then   t    represents 
the  number  of  customers  served  in  the  k —  busy  period  (b.p.) 


and 


they  are  numbered  ^\r--\f    Tk-1  +  ^ '     '  "  '  Tk-1^  '   Next  define 


the  sequence   ^yt,:^  £.  ^    ky 


V1 

Yk  =    I  Wi  ' 

i=T      J 
3    ik-l 

the  sum  of  the  waiting  times  in  the  k —  b.p.   Since  the  queue 

is  stable  for    p  <  1,  we  know  that  W  =»  W  as   n  *  °°  and  wish 

n 

to  estimate   E{W} . 

We  are  interested  in  constructing  controls,  C, ,  which  are 
highly  correlated  with   Z,=Y,-rT,,  but  which  are  not  so  complex 
that  we  can  not  calculate  their  means.   The  controls  we  consider 
are  of  the  form   C.  =  D  -    i./\x,    where  the  r.v.   D,   is  an 
attempt  to  mimic   Y, .   To  compute  the  mean  of   C   we  need,  of 
course,  to  be  able  to  compute    E{t,}.   For  M/G/l  queues  we 
know  that    E{x  }  =  l/(l-p).   For  GI/M/1  queues   E{t  }  =  1/(1-6), 
where  6  is  the  root  inside  the  unit  circle  of   z-E{exp [-u (1-z) u1]}=0, 
where   u,   is  an  interarrival  time.   If  the  queue  in  question  is 
neither  an  M/G/l  or  GI/M/1  queue,  then  the  term   T../y    in   C, 
will  have  to  be  approximated  by  another  r.v.  whose  mean  can 
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be  computed.   The  factor   1/y   multiplying   t,   helps  to  make 
the  variance  reduction  obtained  independent  of  the  scale  parameter 
y.   Also  it  can  be  argued  that  the  term   T,/y   is  then  in  the 
same  units  as   D,  ,  namely,  the  unit  of  time.   Several  alternatives 
for   D    are  listed  below,  indexed  by  a  superscript.   They  are 


(1 


=  xx  =  < 


w0  =  o  , 


wo  +  wi  ' 


Tl =  x 


Ti  1  2; 


D 


(2) 


x1  +  x2  , 


Tl  -  ! 


T1  >  2; 


=  < 


w0  -  0  , 
w0  +  W   , 


.  wQ  +  wx  +  x2  , 


Tl  - 1 
Tl  ■ 2 


ij  >  3; 


and 


D{3)  =  (X+  +  X2)+  =  W2  ; 


D 


(4) 


0  , 


X*  +  (X^  +  x2)+  , 


Tl  =    1 


T,   >_  2 


WQ  =  0   , 

wQ  +  wx  +  w2  , 


Tl  ■ 1 
Tl  ■ 2 
Tl  > 3 
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In  general,  it  is  more  difficult  to  calculate    E{D.   }   as   i 
increases.   On  the  other  hand,  as   i   increases   D     comes 
closer  to   Y,   and  presumably  results  in  more  variance  reduction 
Let   C;[l)  =  D;[l)  -  T-j/u,  i  =  1,2,3. 

In  our  simulation  runs  for  an  M/M/l  queue  the  controls 
C,   ,  C,   ,  and  C,    were  generally  found  to  do  much  better  than 

C,   .   However,  the  more  complicated  controls,  C,     and   C,   , 

(2) 
gave  little  improvement  over   C,   .   Thus  with  both  variance 

(2) 
reduction  and  ease  of  computation  in  mind,  C.     is  our  first 

(2) 
choice.   To  compute   C,    for  the  M/M/l  queue  first  note  that 


P{t1  =  1}  =  P{vQ  <  ux>  =  /   e  Ay  ye~yy  dy  =  (1  +  p)  1 


(ttt7>  f  e~yY  dy        1 

and 


E{xT>  =  E{Xn|Xn  >  0}-  P{X,  >  0}  = 


1     -L-1i-1    ^   -  ^    ^    y (1+p) 


Hence 


E{D1(2)}  =  O  +   P   E{x!  +  X+lxJ  >  0}  =  J*--    [i  +   ,  P  .]  , 
1  1+p     1     2   1         1+p   y    y(l+p) 


(2) 
since   X,   and   X_   are  independent.   The  expectation  of   C, 


is  then  given  by 
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E{C<2)}=  p-lf   P   +   (_^)2  .  jjl 

i       L1+p    i+p    i-pj  * 

Further  discussion  on  the  computation  of  the  mean  of  controls  such 

(2) 
as   C,     for  the  GI/G/1  queue  is  contained  in  the  Appendix.   Com- 
putational results  from  our  simulation  runs  are  contained  in 
Section  3. 


2.3.   Internal  controls  for  Markov  chains 

The  second  example  we  consider  is  a  discrete  time  Markov 
chain.   Assume  now  that  we  are  simulating  an  irreducible,  aperiodic, 
positive  recurrent  Markov  chain,  X  =  {X  :n  ^  0},  with  the  goal 
of  estimating   r  =  E{f(X)}.   The  controls  we  consider  here  are 
of  the  form 

nQA(T1-l) 

(2.6)  C   =     I  [f (X?)  -  r  ]  , 

1      £  =  0  ° 

wtfiere      r„      is    a    guess    of      r,      n        a    fixed   integer,    and 

n   a  (t  -1)   denotes  the  minimum  of   nn   and   (t  -l)  .   Again  our 

motivation  for   C,   is  to  mimic 


Z   =    I      If(X£)  -  r]  =  Yx 


Of  course,  we  still  must  be  able  to  compute   E{C, )   in  order 
to  implement  the  method  of  internal  control  variables.   Let   £ 
be  the  state  space  of   X,  which  we  assume  is  finite,  and   P   the 
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transition  probability  matrix.   Furthermore,  let   g   be  the 

column  vector  with  components   f (i)  -  r    and   .P   be  the  matrix 

0        J 

with  elements 

0  , 


Pkl     I    „         »  -  i  - 


Then  if  we  let   XQ  =  j   and  form  regenerative  cycles  based  on 
the  times  of  return  to  state  j,  it  can  easily  be  shown, using  the 
method  developed  in  Hordijk,  Iglehart,  and  Schassberger  (1976), 

that 

V 
(2.7)  E(C  }  =  (  I         Png)   , 

n=0  J  3 

where   j   is  the  regenerative  state  and  the  subscript  j  means  the 
jth  component  of  the  indicated  vector.   For  two  simple  Markov  chains 
the  theoretical  amount  of  variance  reduction  obtained  using  the 
control   C,   has  been  calculated  for  different  regenerative  states  j 
and  values   n_   and  is  contained  in  the  next  section.   From  (2.6) 
it  is  clear  that  the  larger  the  integer   n_   selected  the  closer 
C-,   comes  to  duplicating   Z,  .   However,  the  larger   n_   the  more 
difficult  is  the  computation  of   E{C, }   contained  in  (2.7). 
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3.   Numerical  Results 

3.1 .   The  M/M/l  queue 

Even  for  the  simple  M/M/l  queue  it  is  in  general  impossible 
to  verify  analytically  what  variance  reduction  will  be  obtained 
via  the  several  internal  controls  listed  in  the  previous  section, 
or  to  get  an  idea  of  the  magnitude  of  the  effect.   For  something 
as  simple  as   C,     it  is  difficult  to  compute  analytically  the 
correlation  between   C,     and   Z,   for  the  M/M/l  queue,  and  this 
is  what  is  required  in  equation  (2.5)  to  find  the  variance 
reduction. 

Thus,  we  resorted  to  simulations  to  verify  the  amount  of 
variance  reduction  obtained  and  the  relative  effectiveness  of 
the  various  controls.   In  the  final  simulations  all  runs  were 
performed  on  an  IBM  System  360/67  computer  using  the  LLRANDOM 
package  (Learmonth  and  Lewis  (1973))  which  generates  random  numbers 
according  to  the  scheme  given  by  Lewis,  Goodman,  and  Miller  (1969) 
and  exponentially  distributed  random  numbers  using  the  Marsaglia 
"rectangle-wedge- tail"  method.   Tests  of  the  random  number 
generator  are  given  in  Learmonth  and  Lewis  (1974)  . 

Of  the  extensive  simulation  runs  performed,  we  give  here 
only  a  summary  of  the  conclusions  and  two  detailed  tabulations 
in  the  case  of  the  most  suitable  control . 

(1)   The  controls   D,   ,  D,     and   D£     do  much  better  generally 
than   D  ' ,  with  little  improvement  over   D,     obtained  by 
use  of   D,     and   D    .   We  say  generally  because  results 
vary  unpredictably  with   X   and   p   and  their  ratio   p. 
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(2)   Subtracting  the  number  of  customers  served  in  a  busy  period 

generally  improves  the  variance  reduction.   By  making  it 

(2) 
dimensionally  stable  as  in   C,     we  obtain  a  "variance 

reduction"  measured  in  terms  of  ratios  of  standard  deviations, 

of  approximately  70%,  uniformly  over   A   and   y.   This  is 

roughly  equivalent  to  halving  the  number  of  cycles  (b.p.'s) 

2 

that  one  must  simulate;  (0.7)   ~  .5.   Much  better  reductions 

can  be  obtained  for  smaller  p  (i.e.  p  =  0.25)  by  specially 

(2) 
designed  controls;  the  point  is  that   C,     works  even  out 

at    p  =  0.99,  where  variance  reduction  is  extremely  important 

Table  1  shows  results  obtained  by  simulating  an  M/M/l 

(2)  .    . 

queue  out  to   n  =  2000  cycles  with  control   C,     and  replicating 

the  simulation  either  m  =  250  or   m  =  100  times  to  estimate  the 

variance  of  the  estimators   r(n),  rCT(n),  and   rCB(n)/  where  we 

drop  the   n   for  convenience.   Here,  we  have  specifically  that 


l   n 
n"  I   Yk 
r^(n)  .  k=lk 


'CB   '     ■>       n 

I 
k=l 


i  I  [xk  +  e(c.[2)  -  e{c{2)})] 


The  estimated  precision  (standard  deviations)  of  the  estimates 
of   E (W)   are  given  in  brackets  under  the  estimates. 

The  results  in  Table  1  are  for   p  =  0.5  (and   m  =  250) 
and  two  values  of   y   and  for   p  =  0.99  (and  m  =  100) 
for  two  values  of   y   commensurate  with  those  given  for 
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p  =  0.5.   The  second,  third  and  fourth  columns  in  the  Table  give 
estimates  of  the  correlations  between  the  control  and  Y,-i,  etc.  from 
which  the  theoretical  variance  reduction  can  be  computed.   They 
are  very  close  to  the  values  given  in  the  next  to  last  column, 
the  observed  variance  reduction,  from  which  we  deduce  that 
estimating   3m   and   3R   affects  the  variance  reduction  only 
slightly.   Overall  there  is  negligible  effect  of  different  values 
of   y   on  the  variance  reduction  obtained  for  the  cases   p  =  0.5 
and   p  =  0.99.   For  the  results   p  =  0.99  given  in  Table  1,  the 
variance  reduction  is  75%,  which  is  about  the  same  as  for 
p  =  0.5.   For  the  case  where  the  control  is  on  the  bottom,  i.e. 
for   t,  the  variance  reduction  is  not  quite  as  good  for  control 
of   Y.   Note  too  that  the  estimated  values  of   E{w}   appear  in 
some  cases  to  be  at  least  three  or  four  standard  deviations  from 
the  true  mean.   This  is  because  the  estimates   r,  r_,m   and   r 

CI  LB 

can  be  seen  from  the  100  replications  to  be  non-normal.   In  other 
words,  for  high   p(0.99),  the  simulation  needs  to  be  taken  out 
further  than  2000  cycles.   In  summary,  the  variance  reduction 
obtained  with  the  internal  control  is  practically  independent 
of  the  scale  factor  y  and  the  traffic  intensity   p. 

3.2.   Two  simple  Markov  chains 

We  turn  now  to  two  simple  Markov  chains  to  give  another 
illustration  of  the  method  of  internal  controls.   The  first 
chain  is  simply  a  recurrent  random  walk  with  state  space 
E  =  {0,1,2,3,4}   and  transition  matrix 
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p  = 


1/2 

1/2 

0 

0 

0 

1/4 

1/2 

1/4 

0 

0 

0 

1/4 

1/2 

1/4 

0 

0 

0 

1/4 

1/2 

1/4 

0 

0 

0 

1/2 

1/2 

Suppose  we  are  simulating  to  estimate   r  =  E{ X} ,  which  can  be 

computed  to  have  the  value  2.   Table  2  contains  the  theoretical 

2       2 

values  of   a  {Z'}/a  {Z,}   and   p(C,,Z,),  where   C,   is  given  by 

(2.6)  for  several  choices  of   n.   and   rn  .   The  entries  in  these 
tables  were  calculated  using  the  methods  in  Hordijk,  Iglehart, 
and  Schassberger  (1976) .   Note  that  the  variance  reduction  is 
largest  when   E{x, }   is  smallest  and   n~   is  largest.   This  also 
yields  the  highest  value  of   p(G,,Z,)   in  keeping  with  our  idea 
of  having   C    mimic   Z,  .   Not  too  much  is  lost  in  variance 
reduction  as  the  guess,   r» ,  departs  from  the  true  value   r  =  2. 
A  variance  reduction  of  more  than  50%  is  attainable,  but  depends 
on  the  state  which  is  used  as  the  regeneration  point. 

The  second  Markov  chain  represents  an   (s,S)   inventory 
model;  see  Crane  and  Iglehart  (1974)  for  further  details.   The 
state  space  £  =  {6, 7,. ..,10}   and  the   P  matrix  is  given  by 
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TABLE  2 

2       2 

Variance  ratio   a   Z'  /a  {Z,}   (top  row) 

and  Correlation   p(C,  ,Z,)   (bottom  row) 


for  the  random  walk  model 


rQ  =  1.5 


1.0 

0.91 

-0.01 

0.30 

0.73 

0.57 

0.52 

0.65 

0.51 

0.40 

0.70 

0.78 

0.91 

0.84 

0.31 

0.40 

0.98 

0.99 

-0.16 

-0.10 

Return  State   i     E{t.}      nn=1       n0=2      nn=5 

0  8       0.97 

-0.16 

1  4        0.82 

0.43 

2  4        0.62 

0.62 

3  4        0.94 

0.25 

4  8        0.97 

-0.16 


r0  "  2 

0  8       0.97  0.98  1.0 

-0.16  -0.12  0.03 

1  4        0.89  0.83  0.70 

0.34  0.42  0.55 

2  4        0.57  0.46  0.35 

0.66  0.74  0.81 

3  4        0.89  0.83  0.70 

0.34  0.42  0.55 

4  8        0.97  0.98  1.0 

-0.16  -0.12  0.03 
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TABLE    3 

2  2 

Variance    ratio      a    {Z'}/a    {  Z,  }        (top   row) 

and  Correlation   p(C, ,Z.)   (bottom  row) 

for  (s,S)  Inventory  Model 


ro  =  8 
Return  State  i         E   «     nn=1         nn=2       nn=4 


8 

9 

10 


5.55 

0.75 

0.65 

0  .49 

0.51 

0.59 

0.72 

5.72 

0  .82 

0.67 

0  .49 

0.43 

0.57 

0  .72 

6.27 

0.89 

0.76 

0  .57 

0.33 

0.49 

0.66 

7.21 

0.95 

0.86 

0.69 

0.21 

0.38 

0  .56 

2.88 

0.95 

0  .81 

0  .42 

0.22 

0.44 

0.77 

r0    "    9 

5.55 

0  .75 

0.75 

0  .78 

0.51 

0.50 

0.47 

5.72 

0.88 

0.80 

0  .73 

0.34 

0.45 

0.52 

6  .27 

0.93 

0  .83 

0.71 

0.27 

0.41 

0.54 

7.21 

0.95 

0.86 

0.72 

0.23 

0.37 

0.53 

2.88 

0.79 

0.57 

0.24 

0.46 

0.66 

0.87 
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Again  we  wish  to  estimate   E{X}  ■  8.297.   Table  3  contains  the 
results  analogous  to  those  for  the  random  walk  (Table  2) .   The 
same  general  observations  about   E{t,}   and   n_   made  for 
Table  2  seem  to  hold  here  as  well.   Note  however  that  the 
results  are  sensitive  to  the  choice  of   rQ .   Again  a  50% 
reduction  in  variance  can  be  attained. 


4 .   Conclusions  and  Extensions 

We  have  been  able  to  obtain  a  50%  variance  reduction  using 
internal  control  variables,  for  the  regenerative  estimate  of  the 
limiting  value  of  the  mean  waiting  time  in  an  M/M/l  queue.   This 
reduction  is  obtained  uniformly  over  all  parameter  values.   It  is 
fairly  certain  that  the  technique  will  work  well  with  any  GI/G/1 
queue. 

Internal  control  variables  can  be  easily  used  with 
discrete  time  Markov  chains.   The  examples  used  in  this  paper 
showed  that  a  variance  reduction  of  50%  is  attainable.   This  figure 
is  likely  to  vary  widely  with  the  particular  Markov  chain. 
Continuous  time  Markov  chains  and  semi-Markov  processes  can  be 
handled  in  the  same  way  using  the  discrete  time  method  of 
Ilordijk,  Iglehart  and  Schassberger  (1976)  . 

Another  method  of  internal  stratified  sampling  was  also 
investigated.   This  method  produced  little  overall  variance 
reduction  despite  considerable  effort. 
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Finally  we  note  that  it  is  possible  to  apply  the  idea  of 
internal  controls  to  the  classical  sample  average  estimate  over 
a  realization  of  fixed  length   m.   Thus  in  estimating   E(W) 
in  the  GI/G/1  queue  we  have 


m-1 

W(m)  =  -   I   W.  , 
j-0 


(4.1) 


which  may  be  written  as 

-,   N(m) 

W(m)  =  -  {   T   Y,  +  Y'  ,  w ,}  ,  (4.2) 

m  L  ,^,   k     n  (m)  +1J 

where   N(m)   is  the  number  of  completed  busy  periods  in  the  queue 

in   [0,m],  Y,   as  before  is  the  sum  of  the  waiting  times  in  the 

kth  cycle,  and   Y ' ,  >  .   is  the  sum  of  the  waiting  times  in  the 

N  (m)  +1  J 

last,  incomplete  cycle. 

A  central  limit  theorem  for  W(m)   holds  for  which  the  variance 

2 
term  is  proportional  to   E{Z,  },  and  which  is  now  estimated 

from  the  random  number   M (m)   of  cycles.   Control  is  applied 
to  the   Y,  ' s  in  (4.2)  just  as  it  is  applied  in  the  ratio 
estimator   r__,(n).   Call  this  estimator   W  (m)  .   For  the  M/M/l 
queue  the  variance  reduction  observed  in  the  simulations  was 
the  same  for  all  controls   C,     with  the  ratio  estimator  and 
with   W  (m) . 

The  main  reason  for  considering  W  (m)   is  that,  while 
it  loses  the  advantage  of  being  an  estimator  using  a  fixed  number 
of  i.i.d.  random  variables,  one  can  apply  the  classical  external 
controls  to   W  (m) .   Thus  one  could  use  the  difference  of  the  sum 
of  the   m   arrival  times  and  the   m   service  times  to  control 
W  (m) .   We  have  not  yet  tried  this  idea. 
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APPENDIX 

To  implement  the  internal  control  techniques  for  a 
GI/G/1  queue,  certain  theoretical  parameters  are  required.   In 
this  appendix  we  shall  indicate  the  values  of  these  parameters 
in  so  far  as  they  can  be  calculated.   These  values  are  either 
well-known  or  easily  calculated.   For  a  reference  to  the  known 
formulas  see  Cohen  (1969). 

We  begin  with   E{x  },  the  expected  number  of  customers 

served  in  a  busy  period.   For  the  general  GI/G/1  queue  recall 

that  we  let   X   =  v   ,  -  u    and   S   =  X,  +  • • •  +  X  ,  for   n  >  1, 
n    n-1     n        n     1  n         — 

with   S„  =  0 .   Then    t,  =  inf{n  >  0:S   <  0}.   The  general 
0  1  n  —  3 

expression  for   E{t,}   is  given  by 


E{t,}  =  exp{  V   n  1P[S   >  0]} 
1         u  ->  n 

n=l 


an  impossible  expression  to  evaluate  in  general.   Another  useful 


expression  for   E{x,}   is 


(A.l)  E{Tl)  =  1/P{W  =  0}  , 


where   W   is  the  stationary  waiting  time.   In  the  special  case 
of  the  M/G/l  queue,  however,  we  have 


E{Tn}  =  (1  -  p)  1 
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Now  for  the  queue  GI/M/1  we  can  use  (A.l)  and  the  stationary 
distribution  of  the  embedded  Markov  chain  to  conclude  that 


E{t1>  =  (1  -  6)  X 


where   6   is  the  root  inside  the  unit  circle  of 

z  -  U{y  (1-z) }  =  0 

-su1 
with   U(s)  =  E{e     },  Re  s  >_  0 ,  and  where   v0   is  an  exponential  (y) 

r.v.   It  is  easy  to  check  that   6  =  p   for  M/M/l  queues.   Daley 

(1975)  has  recently  proposed  the  approximation  to   6   given  by 


6  =  a,(l-p)2  +  2(l-b  1)p  +  (2b  1-l)p2/ 


2 

where   a,  =  P{u  =0},  E{u,  }  =  1,  and   b  =  E{u-,}.   This  approximation 

gives  good  results  in  a  number  of  examples  calculated  by  Daley 
(1975)  and  may  be  useful  for  the  purposes  of  this  paper. 

Next  we  turn  to  the  computation  of   P{t,=1)   and   P{x  =2} 
For  the  GI/G/1  case  we  have 


p{Tl=i}  =  P{S1  <_   0} 


and 


P{t1  =  2}  =  P{S1  >  0,  S2  <_  0}, 


25 


both  of  which  can  be  worked  out  with  a  little  effort.   For  the 
M/M/l  queue 


and 


P{x1=l)  =  (1+p)'1 


P{x1=2>  =  p(l+p)"3  . 


For  the  M/G/l  queue 

P{t1=1}  =  V(A)  , 


"Avo 

where   V(A)  =  E(e     ),  while  for  the  GI/M/I  queue 


P{t1=1)  =  1  =  U(y)  , 


whe 


re   U(s)   is  given  above.   For  the  M/E./l  and  E./M/l  queues 


the  value   P{t,=2}  can  be  calculated  with  some  effort.   As  these 
expressions  are  cumbersome  they  shall  be  omitted. 

Finally  we  give  various  partial  expectations  which  are 

(4) 
needed  for  computing  the  means  of  internal  control   C,   .   Namely, 


E{S*  +  S*/  tx  >_  2}  =  EfS^  S1    >  0}  +  EfS*,  S,  >  0} 
and 

E{t1'  t1  -  2}  =  e{ti}  ~  Ht.^1}. 


Here  the  symbol   E{X,A}  =  E{Xl  } ,  where   X   is  a  r.v.,  A   an 
event,  and   1    the  indicator  function  of   A.   In  the  special 
case  of  the  M/M/l  queue, 
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E{S1#    S    >    0}    =    p/{y (1+p) }     , 


■<•  si>°»  =[2(r^)2  +  T^p]^  - 


and 


E 


{xl'    ti    1  2}    =    2P/{  (1"P)  d  +  P>>     • 
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